#########################################################
# replication Rasmussen 2010 [origin]
#########################################################
library(readstata13)

# downloaded from http://journals.sagepub.com/doi/suppl/10.1177/1465116510388675
dat_r = read.dta13("replication_Rasmussen/EUP388675-supplementary_material.dta")
dat_r$var15 = as.numeric(as.character(dat_r$var15))
dat_r$var5 = as.integer(dat_r$var5)
dat_r$var5[which(dat_r$var5==1)] = 0
dat_r$var5[which(dat_r$var5==2)] = 1

dat_r$var16 = as.integer(dat_r$var16)
dat_r$var16[which(dat_r$var16==1)] = 0
dat_r$var16[which(dat_r$var16==2)] = 1

dat_r$var10 = as.integer(dat_r$var10)
dat_r$var10 = dat_r$var10 - 1 

dat_r$var15 = dat_r$var15 - 1998

dat_r1 = subset(dat_r,dat_r$var13==0)

# Covariates: 
# 
# Policy distance rapporteur-EP median
# Big party group
# Country coherence of key negotiators
# Closeness of working relationship (calendar year)
# Urgency of legislation: existing legislation expires; No existing legislation
# Workload (pending legislation during first reading)
# Issue areas consulted (no. of consultative committees)
# New act
# Directive
# Inter-institutional disagreement
# 

# Focusing on risk of delegation

dat_r1 = subset(dat_r1,select = c(var16,var6,var3, var5,var10,var15,var1,var4,var11,var14,
                                  var7,var17))


dat_r1 = dat_r1[complete.cases(dat_r1),]

# Replicate findings in Model III of Table 2
model1 = glm(var16 ~ var6 +var3+ var5 + var10+ var15+ factor(var1)+
               var4+ var11+ var14+ var7+ var17,
             data=dat_r1,family = binomial(link="logit"))
summary(model1)
pred = predict(model1,dat_r1,type="response")
pred_y = ifelse(pred>0.5,1,0)
pred_ratio = sum(pred_y == dat_r1$var16)/385 # 76.6%

###############################
# counterfactuals

names(dat_r1)
head(dat_r1,40)
dat_r2 <- with(dat_r1, data.frame(var6 = seq(0,2,length.out = 1000), 
                                  var3=seq(0,2,length.out = 1000), 
                                  var5=1,
                                  var10=median(var10), 
                                  var15=median(var15), 
                                  var1="Yes", 
                                  var4=median(var4,na.rm=T), 
                                  var11=median(var11, na.rm=T), 
                                  var14="Yes", 
                                  var7="Yes", 
                                  var17="Yes"
                                  )
               )

preds = predict(model1,dat_r2,type="response", se.fit=TRUE)

predm = preds$fit
predl = preds$fit - (1.96*preds$se.fit) # lower bounds
predu = preds$fit + (1.96*preds$se.fit) # upper bounds







